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NOMENCLATURE 


Dynamic  viscosity  coefficient 

Velocity  of  sound 

Friction  coefficient 

Transformed  coordinate 

Body  diameter 

Velocity  functions 

Parameter  defined  in  Eq.  (2.8) 

Convective  heat  transfer  coefficient 

Parameters  defined  in  Eq.  (2.17) 

Thermal  conductivity 

Mach  number 

Nusselt  number 

Pressure 

Prandtl  number 

Heat  flux 

Reynolds  number 

Parameters  defined  in  Eq.  (2.7) 
Relative  enthalpy  difference 
Temperature  functions 
Temperature 
Tangential  velocity 
Coefficients  of  velocity  series 
Boundary  layer  edge  velocity 
Coordinates  along  body  surface 


Transformed  coordinate 


Coordinate  normal  ro  body  surface 
Angle 

Velocity  parameters  defined  in  Eqs.  (2 

Specific  heat  ratio 

Similar  independent  variable 

Dynamic  viscosity 

Kinematic  viscosity 

Density 

Shear  stress 


Subscripts 


Stagnation 
Boundary  layer  edge 


.7)  and  (2.9) 
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I.  INTRODUCTION 

Aerodynamic  heating  of  a blunt  body  la  a result  of  flow  of  air  at 
high  speed  about  it.  Direct  compression  and  Internal  friction  at  and 
near  the  stagnation  regions  of  forward  surfaces  of  the  body  convert  the 
kinetic  energy  of  motion  into  heat  within  the  boundary  layer  of  air 
which  surrounds  the  body.  Consequently,  aerodynamic  heating  probleno 
dictate  the  design  of  blunt  bodies,  not  only  for  structural  reasons, 
but  also  because  of  thermal  problems  associated  with  protecting  vital 
internal  components,  such  as  electronic  packages  which  are  usually 
located  in  the  forward  part  of  the  body. 

Previous  investigations  concerning  aerodynamic  heating  problems 
have  concentrated  in  the  area  of  high  supersonic  or  hypersonic  flows 
at  high  altitude  low  density  atmospheric  flight,  such  as  in  entry  or 

r ,* 

re-entry  cases. I 1]  However,  current  interest  in  blunt  bodies  which 

fly  at  moderate  supersonic  speeds,  but  in  the  dense  low  altitude  atmos- 
phere, call  for  further  investigation  in  aerodynamic  heating  problems 
characterized  by  relatively  lower  Mach  and  Reynolds  numbers . The 
objective  of  this  study  is  primarily  to  obtain  heat  transfer  coefficients 
and  recovery  factors  along  the  forward  surfaces  of  the  body  to  allow 
for  coupling  these  parameters  into  a complete  heat-transfer  analysis 
Inside  the  body.  This  information  can  be  obtained  by  solving  s three- 
dimensional,  compressible  boundary  layer  around  a body  with  a blunt 
nose,  which  may  or  may  not  have  an  additional  rotating  speed  complication. 

_ 

Bracketed  numbers  refer  tc  entries  in  REFERENCES. 
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The  basic  nonlinear  partial  differential  aquations  [2]  which 
govern  the  no f—u*.  of  a steady,  axisyametric,  compressible  nonrotating 
1 sains r boundary  layer  flow  about  a body  of  r rvolution  >tave  been 
transforaed  into  a aore  convenient  fora  by  a nod if led  Illingworth- 
Stewsttson  transforastion  [3].  A special  procedure  to  relate  the 
physical  sensible  external  flow  conditions  with  the  transforaed  ones 
was  presented.  Siallarity  variables  were  found  by  applying  the  system- 
atic one-paraaeter  group  theory.  The  siapllfled  governing  aquations 
were  then  transforaed  again  by  the  well-known  siailar  analysis  [3]. 

A perturbation  scheae  based  on  the  transforaed  coordinates  was  con- 
structed to  render  a series  of  coupled  nonlinear  ordinary  differential 
equations  which  are  readily  solved  by  standard  numerical  Integration 
subroutines  to  provide  the  desirable  flow  property  distributions,  in- 
cluding heat  transfer  rates.  The  purposes  of  this  report  are  to  present 
solutions  to  the  differential  equations  and  to  examine  velocity  and 
teaperature  profiles  within  the  boundary  layers  as  well  us  shear 
stresses  and  hest  transfer  rates  at  the  surface  of  the  body. 


2.  ANALYSIS 


2 . 1 Governing  equations 

The  basic  nonlinear  partial  differential  equations  which  describe 
mass,  momentum  and  energy  transport  for  steady,  axisymmetric,  compress- 
ible, nonrotating  laminar  boundary  layer  flow  about  a body  of  revolution 
have  been  transformed  into  a series  of  coupled  nonlinear  ordinary  dif- 
ferential equations.  Details  describing  this  analysis  are  presented 
elsewhere  [3]  and  are  not  cited  here.  For  convenience,  however,  only 
the  ordinary  differential  equations  are  presented.  These  equations  are 
given,  after  some  modification,  from  those  previously  reported  [3] 
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with  corresponding  boundary  conditions 
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n = 0 : 


fl  " 0 g3  “ 0 g5  ■ 0 h5  - 0 


f [ - 0 «3  - 0 85  - 0 h5  - 0 


wj  - 0 


z.  **  0 
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Solution  of  these  differential  equation?  with  associated  boundary  con- 


ditions is  presented  in  Chapter  3 with  results  discussed  in  Chapter  4. 
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2.2  Boundary  layer  characteristics 

Once  solutions  to  the  ordinary  differential  equations  given  in 
Sec.  2.1  are  available,  quantities  such  as  velocity  and  temperature 
profiler,  within  the  boundary  layer  as  well  as  shear  stresses  and  heat 
transfer  rates  at  the  surface  of  the  body  can  be  evaluated.  The  purpose 
of  this  section  is  to  present  expressions  which  may  be  utilized  to 
evaluate  these  quantities  as  well  as  others  that  are  generally  employed 
to  describe  boundary  layer  characteristics.  It  is  convenient  to  intro- 
duce certain  dimensionless  variable  which  allow  presentation  of  results 
and  discussion  to  be  considerably  simplified.  Examination  of  definitions 
and  expressions  for  transformed  variables,  stream  function  and  relative 
enthalpy  difference  as  presented  in  Ref.  [3]  reveals  that  the  following 
dimensionless  variables  can  be  defined  in  terms  of  the  previously  intro- 
duced parameters. 


x " x/D 


d » x/bD 


u3bV 


i + 3) 
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4 4 
u,b  D 
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i 2(28g2  + 72g  + 45) 
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where 
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0 is  expressed  as 


$ - M 


[”a[(Y  - i)  m*  + 2) 

, , v + 1 - » M-  + 2 

| (y  + i) 

1 + 2 2 

2Y  M‘  - (y  - 1) 

-U/<Y-1>] 


(2*9) 

Employing  the  above  dimensionless  variables,  the  distance  along  the 
surface  measured  from  the  tip  cf  the  blunt  body  is  expressed  as 


x - d + 3^  4^  + g (2.10a) 

m 

For  positions  near  the  tip  of  the  body,  d is  small  and  Eq.  (2.4)  reduces 
to 


x ■ cl  (2.10b) 

Because  of  the  method  employed  in  the  analysis  to  describe  the  boundary 
layer  edge  velocity,  x is  limited  to  values  less  than  about  0.611  which 
corresponds  to  an  angle  a as  defined  in  Ref.  [3]  of  35°.  The  maximum 
value  of  the  transformed  coordinate  d is  dependent  on  the  limiting 
value  of  x.  In  terms  of  x,  d can  be  written  as 

d-x-x^  4^  + x5J2g  (2.11) 

In  expressing  Eqs.  (2.10)  and  (2.11),. ud  to  ne  fifth  order  results  were 
retained  in  the  perturbation  scheme  [3].  Similar  expressions  for  the 
distance  measured  normal  to  the  body  surface  can  be  developed  but  are 
not  presented  here. 
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In  most  applications,  the  tangential  velocity  component  within  the 
boundary  layer  is  of  interest  since  it's  gradient  determines  the  drag 
experienced  by  the  body  as  it  passes  through  a gas.  Consequently,  only 
results  for  the  tangential  velocity  are  presented.  Expression  for  the 
normal  velocity  component  may  be  developed  using  the  definitions  as 
reported  in  Ref.  [3].  Utilizing  the  definitions  for  the  stream  function 
as  well  as  for  the  boundary  layer  edge  velocity,  the  normal  velocity 
component  is  expressed  as 


fl  + 81  83  + ^<*2  85  + 8l  hp 

1 + d^  s^  + 3^  Sj 


(2.12a) 


where  represents  the  boundary  layer  edge  velocity.  Velocity  functions 
f1#  g-j.  g5  and  h^  are  solutions  to  the  ordinary  differential  equations 
and  are  understood  to  be  dependent  on  the  similar  independent  variable 
r\.  At  a sufficiently  large  value  of  rj  corresponding  to  the  boundary 
layer  edge,  the  velocity  ratio  in  Eq.  (2.12a)  attains  a value  of  unity. 
For  small  values  of  d,  this  ratio  is  given  by 


“ - f£(n)  (2.12b) 

which  is  a solution  to  Eq.  (2.1).  At  the  surface  of  the  body  where 
n ■ 0,  the  velocity  ratio  is  zero  as  can  be  observed  by  imposing  the 
boundary  conditions  given  in  Eq.  (2.5). 

Temperature  profile  within  the  boundary  layer  is  determined  from 
the  definition  of  the  relative  enthalpy  difference  [3].  After  some 


manipulations  and  recognizing  the  dimensionless  variables,  the  gas 
temperature  normalized  with  the  free  stream  stagnation  temperature  is 
written  as 
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sr*  - (S  + 1)  - <3x)2 

r0  1 U1 


(2.13a) 


The  relative  enthalpy  difference  S Introduced  in  Eq.  (2.13)  is  given 
in  terms  of  solutions  to  the  differential  equations  as 


S - SQ  + d2  sL  z2  + d4(s2  z^  + s2  w £ 


(2.14) 


where  temperature  functions  Sq,  z^,  z^  and  w^  are  functions  of  n*  On 
the  surfer*  of  the  bcdy,  Eq.  (2.13)  reduces  to 

T 

f2  - • + 1 «-15) 

iQ  W 

where  Tw  is  surface  temperature  which  may  vary  with  x.  Values  of  Sw 
less  than  zero  correspond  to  cooled  surface  ?w/Tq  < 1 end  greater  than 
zero  to  a heated  surface  T^/Tq  > 1.  At  the  boundary  layer  edge,  S ■ 0 
by  boundary  conditions  cited  in  Eq.  (2.6)  and  the  temperature  is 


(2.16) 


where  subscript  "1"  refers  to  the  boundary  layer  edge.  Thus,  the 
boundary  layer  edge  temperature  is  just  a function  of  x and  free  stream 
Mach  number.  In  the  analysis  [3],  the  Prandtl  number  was  assigned  a 
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value  of  unity.  Hence,  the  adiabatic  wall  boundary  condition  yields 
S - 0 within  the  boundary  layer  which  implies  that  the  temperature 
distribution  is  solely  dependent  on  the  tangential  velocity  distribution 
for  a given  x position  [4].  This  can  be  observed  by  substituting  a 
value  of  zero  for  S in  Eq.  (2.13a).  In  addition,  the  wall  tenperature 
for  Pr  ■ 1 with  adiabatic  boundary  condition  eq>ials  the  free  stream 
stagnation  temperature  and  the  recovery  factor  is  unity  [5].  For 
sufficiently  small  values  of  cl,  Eq.  (2.13)  may  be  written  as 

- S.(n)  + 1 (2.13b) 

l0  u 


where  it  vac  further  assumed  that  the  velocity  is  small.  This  is  justi- 
fiable  since  small  values  of  d correspond  to  the  stagnation  region  of 
the  blunt  body  where  the  velocities  are  small. 

Shear  stress  at  the  body  surface  indicates  the  drag  on  the  body 
and  is  evaluated  from  the  following  expression 
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3u 
^w  3y 
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“VbwoYpoB 


Bx 


M®  Ifr  * **  81  83'  + 85  * B1  h5,)]n  ■ 0 

To)  1 + d2  s1  + d4  7? 


(2.17a) 


Second  derivatives  of  the  velocity  functions  are  thus  related  to  the 
shear  stress.  Equation  (2.17)  reduces  to  the  following  for  small  values 
of  d 


T 


W 


^by0YP0B  Bx  fj'(0) 


(2.17b) 
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.Another  quantity  of  Interest  Is  skin  friction  coefficient  given  by 


r 2 (S  + 1)  v7‘ 

w N w ' 0 

' 3 V“x°2 


tfr  ^ S1  8j'  + d*<»2  >5'  + «1  ■ 0 

-2  — -2  — 2 

(1  v cT  sx  + <T  s2) 


(2.18a) 


For  small  values  of  d,  skin  friction  coefficient  is 


2(s«+1)  ST  f.. 

i Vux  °2  1 


(0) 


(2.18b) 


Finally,  the  Reynolds  number-skin  friction  parameter  is  expressed  as 


Cf  V^7  , JTjuTT 

2 if  d In  x 


f£l’  + 8i  83  + ^ g5  + S1  h5^n  - 0 
n „ A _ 3/2 

(1  + d + 3 s2) 


where 


(2.19a) 


d In  d _ 1 + d2  1 g/3  + d4  i?  k(7k  + 3)/30 

d In  x , *2t  . "4  —2  ...  , ,, 

1 + d j g + d J g(/g  + 3)/6 


(2.20) 
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The  Refolds  number  Introduced  In  Eq.  (2.19)  la  defined  as 

Re  - U,  x/v  (2.21) 

where  all  properties  are  evaluated  at  the  body  surface.  Hear  the 
stagnation  region  of  the  blunt  body,  this  parameter  assumes  the  form 

» f['(0)  (2.19b) 

Heat  flux  at  the  body  surface  is  evaluated  from 


q — k 


ar 


w 3y  j 


s’(0) 


(2.22a) 


where  S'(0)  is  found  by  taking  derivative  of  Eq.  (2.14)  with  respect  to 
n and  evaluating  at  n * 0.  Thus,  heat  transfer  is  related  to  fi*.8t 
derivatives  of  the  temperature  functions.  For  small  d values,  Eq.  (2.22) 
is  written  as 


(2.22b) 


Positive  and  negative  values  of  heat  flux  imply  heating  and  cooling  of 
the  surface,  respectively.  The  boundary  condition  with  S'(0)  ■ 0 yields 
a zero  heat  flux  which  corresponds  to  adiabatic  or  insulated  wall.  A 
convective  heat  transfer  coefficient  can  be  defined  as  follows 


q - h(Tw  - T0) 


(2.23) 
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where  the  free  stream  stagnation  temperature  has  been  employed.  Several 
other  definitions  [4,5]  employ  the  adiabatic  wall  temperature  which  is 
the  temperature  acquired  by  an  adiabatic  surface.  However,  for  the 
existing  analysis  with  Pr  - 1,  the  adiabatic  and  free  stream  temperatures 
are  identical.  The  local  Nusselt  number  can  then  be  written  as 


„ hx 
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V Jd  )n  ? 
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-S1CQ>. 
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(2.24a) 


This  expression  reduces  to  the  following  form  for  small  values  of  d 
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(2.24b) 


The  dependency  of  the  Nusselt  number  on  the  Reynolds  number  as  displayed 
in  Eq.  (2.24)  is  similar  to  that  observed  for  a flat  plate,  cylinder  or 
sphere  for  laminar  boundary  layer  flow  [6]. 

A final  parameter  that  is  useful  in  examining  boundary  layer 
characteristics  is  the  Reynolds  analogy  parameter  which  for  the  present 
analysis  acquires  the  form 
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For  small  vaues  of  d,  this  parameter  reduces  to 
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! I 


The  Reynolds  analogy  parameter  illustrates  the  interrelationship 
between  fluid  friction  and  heat  transfer  processes.  If  the  Prandtl 
number  is  included  in  this  parameter,  then  the  dimeasionlee*  grouping 
of  Nu/Rc^Pr  is  known  as  the  Stanton  number.  Laminar  boundary  leyer  on 
a flat  place  with  zero  pressure  gradient  yields  a value  of  one  for  the 
Reynolds  analogy  parameter  [7]. 


14 


3.  METHOD  OF  SOLUTION 

The  method  employed  U"  obtain  solutions  to  the  ordinary  differen- 
tlal  equations  as  given  in  Sec.  2.1  was  a fourth-orciar  Kunge-Kutta 
integration  scheme  using  double  precision  arithmetic  on  IBH  360/65 
digital  computer  system.  It  was  found  advantageous  to  first  solve  the 
set  of  equations  given  in  Eq.  (2.1)  with  appropriate  boundary  conditions. 
Using  these  results,  solutions  to  Eqs.  (2.2)  and  (2.3)  were  acquired. 
Finally,  the  results  for  Eqs.  (2.1)  and  (2.2)  were  employed  to  obtain 
solutions  to  Eq.  (2.4).  i'he  technique  for  solving  each  set  of  equations 
is  now  outlined  with  additional  information  for  solving  ordinary  dif- 
ferential equations  of  the  boundary  layer  type  available  elsewhere  [8,9]. 
Since  the  integration  scheme  requires  knowledge  of  values  for  the  func- 
tions as  well  as  their  derivatives  at  n * 0,  initial  guesses  for  the 
unknown  derivatives  (for  example,  in  Eq.  (2.1),  f'^O)  and  Sq(0)  are 
unknown)  were  made  and  the  integration  carried  out  to  some  value 

(initially  2)  where  boundary  conditions  specified  in  Eq.  (2.6)  must  be 
satisfied.  If  these  boundary  conditions  were  not  met,  then  the  guesses 
for  the  derivatives  must  be  adjusted  and  the  integration  repeated.  The 
method  employed  to  obtain  new  estimates  for  the  derivatives  is  attri- 
buted to  Nachtschelm  and  Swlgert  [9].  Upon  satisfaction  of  the  boundary 

conditions  for  the  particular  value  of  n , it  was  then  necessary  to 

max 

establish  if  n x corresponded  to  a sufficiently  large  value  as 

required  by  Eq.  (2.6).  n was  then  increased  (for  example,  next 

max 

value  was  4)  and  the  integration  scheme  as  well  as  adjustment  of  values 
for  the  unknown  derivatives  repeated  until  the  boundary  conditions  were 
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again  satisfied.  This  procadura  vaa  continued  until  values  for  the 

unknown  derivatives  at  n ■ 0 did  not  vary.  In  all  cases,  the  ssxlnis 

value  for  n was  6 where  all  boundcry  conditions  were  generally  within 
-12 

10  of  the  required  values.  An  integration  step  size  of  0.005  was 
found  to  give  sufficiently  accurate  results  of  at  least  eight  decimal 
digits  for  all  initial  darivatlvaa  and  reasonable  computational  tines. 

A listing  of  the  digital  computer  program  to  obtain  valuea  of 
unknown  initial  derlvativas  is  given  in  Appendix  A.  Also,  a program 
which  uses  these  results  to  generate  the  values  for  all  functiona  at 
different  n values  for  listing,  plotting  and  analysis  purposes  Is 
supplied.  Results  from  this  numerical  method  are  presented  in 
Chapter  4. 
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4.  DISCUSSION  OF  RESULTS 

4.1  Solution*  of  governing  equations 

Utilising  the  method  of  solution  ac  discussed  in  Chapter  3, 
solutions  for  the  unknown  vsluss  of  derivative*  st  ths  body  surface 
ware  acquired  for  several  values  of  vail  enthalpy  difference  S from 
-1  to  6 and  results  presented  in  Table  4.1.  Employing  these  values  as 
well  an  those  in  Eq.  (2.5)  as  ini tie?  conditions  for  the  Runge-Kutta 
integration  scheme,  velocity  and  temperature  functions  were  evaluated 
for  different  values  of  the  similar  independent  variable  and  are  pre- 
sented in  Figs.  4.1,  4.2,  4.3  and  4.4  for  Eqs.  (2.1),  (2.2),  (2.3)  and 
(2.4),  respectively.  Only  first  and  second  derivative  results  for  the 
velocity  functions  are  displayed  since  these  correspond  to  velocity 
and  shear  stress,  respectively.  dentified  in  Sec.  2.2,  the  various 

quantities  used  to  describe  boundary  layer  characteristics  are  expressed 
in  terms  of  the  functions  f ^ and  Sq  for  positions  near  the  stagnation 
point.  These  functions  are  shown  in  Fig.  4.1.  Only  results  are 
illustrated  for  values  of  n up  to  4 where  it  can  be  observed  in  the 
graphs  that  the  boundary  conditions  specified  by  Eq.  (2.6)  are  adequately 
satisfied.  Several  general  comments  can  be  made  concerning  results 
presented  in  these  figures.  First,  tangential  velocities  within  the 
boundary  layer  may  exceed  the  boundary  layer  edge  velocity.  This  can 
be  observed  by  recognising  the  greater  than  unity  values  shown  by  the 
first  derivatives  of  the  velocity  functions.  Furthermore,  the  second 
derivatives  of  the  velocity  functions  increase  as  wall  enthalpy  dif- 
ference increases,  and  the  shear  stress  is  expected  to  exhibit 
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a similar  trend.  The  adiabatic  boundary  condition  yields  S ■ 0 
throughout  the  boundary  layer.  This  Implies  that  the  total  energy 
(Cp  Tq)  remains  constant  within  the  boundary  layer  and  the  surface 
is  at  the  free  stream  stagnation  temperature.  Finally,  heating  of  the 
surface  by  the  gas  occurs  for  wall  enthalpy  differences  less  than  zero 
and  cooling  for  wall  enthalpy  differences  greater  than  zero. 

4.2  Results  for  boundary  layer  characteristics 

Velocity  and  temperature  profiles  within  the  boundary  layer  for 
several  locations  along  the  body  and  for  various  wall  enthalpy  dif- 
ferences are  illustrated  in  Fig.  4.5  for  a Mach  number  of  1.5  and 
specific  heat  ratio  of  1.4.  Results  for  x * 0 correspond  to  the  stag- 
nation region  where  the  boundary  layer  edge  velocity  is  zero.  At  this 
location  there  wculd  be  no  hydrodynamic  boundary  layer.  Positions 
slightly  removed  from  the  stagnation  point  have  velocity  profiles 
represented  by  those  for  x - 0 as  Illustrated  by  Eq.  (2.12b).  Tangen- 
tial velocities  greater  than  the  boundary  layer  velocities  are  dis- 
played for  the  higher  values  of  wall  enthalpy  difference.  Velocity 
ratios  greater  than  unity  are  attributed  to  the  Increase  in  volume 
imparted  to  the  gas  due  to  the  high  wall  temperatures.  The  gas  of 
lower  density  is  accelerated  by  the  pressure  in  spite  of  it  being 
decelerated  by  viscous  forces.  A thermal  boundary  exists  near  the 
stagnation  region  for  wall  temperatures  different  from  the  free  stream 
stagnation  temperature  and  grows  along  the  body.  The  decrease  in  tem- 
perature ratio  along  the  body  is  a result  of  an  increase  in  the  amount 
of  stagnation  energy  being  transformed  into  kinetic  energy.  The 
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boundary  layer  edge  velocity  Increases  from  the  stagnation  point  with 
a corresponding  reduction  in  boundary  layer  edge  temperature.  Results 
for  nonunifom  wall  t miperatures  can  be  obtained  from  those  presented 
in  Fig.  4.5  by  specifying  the  wall  temperature  for  each  position  and 
selecting  the  corresponding  curve.  Boundary  layer  edge  velocity  and 
temperature  results  are  not  influenced  by  the  wall  temperature. 

As  previously  mentioned,  surface  shear  stress  is  related  to  the 
drag  experienced  by  a body  as  it  passes  through  a gas.  Representative 
values  for  shear  stress  normalized  with  respect  to  the  factor  of 
/bUgY Pq3  are  illustrated  in  Fig.  4.6  as  a function  of  distance  along 
the  body  for  a free  stream  Mach  number  of  1.5  and  specific  heat  ratio 
of  1.4.  The  maximum  distance  along  the  body  for  which  the  analysis  [3] 
applies  is  limited  by  the  applicability  of  the  method  to  describe  the 
boundary  layer  edge  velocity.  Results  for  a particular  value  of  wall 
enthalpy  difference  correspond  to  an  isothermal  surface.  Fo*  noniso- 
thermal  surfaces,  similar  results  can  be  obtained  from  those  presented 
in  Fig.  4.6  provided  the  distribution  of  wall  enthalpy  difference  along 
the  body  is  specified.  Results  illustrate  that  wall  shear  stress  is 
lower  when  the  surface  is  cooled.  This  is  attributed  to  wall  dynamic 
viscosity  (y  ^ T ) and  second  derivatives  of  velocity  functions  (see 
Table  4)  exhibiting  smaller  values  on  a cooled  surface.  Near  the 
stagnation  region,  shear  stress  increases  almost  linearly  with  distance 
as  also  can  be  observed  from  Eq.  (2.17b).  The  shear  stress  is  found 
not  to  be  a strong  function  of  wall  temperature.  For  example  a twofold 
increase  of  wall  enthalpy  difference  from  1.0  to  2.0  yields  only  about  a 
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20X  increase  in  shear  stress  at  a location  of  0.4.  The  decrease  in 
shear  stress  with  increasing  distance  is  believed  to  be  a result  of  the 
boundary  layer  beginning  to  separate  from  the  body.  However,  further 
investigation  is  needed  to  establish  the  validity  of  this  conjecture. 
Results  for  the  Reynolds  number-skin  friction  parameter  as  defined  in 
Eq.  (2.19a)  are  displayed  in  Fig.  4.7  as  a function  of  distance  for  a 
Mach  number  of  1.5  and  specific  heat  ratio  of  1.4.  Values  of  this 
parameter  for  points  near  the  stagnation  region  are  given  by  values 
of  f^'(0)  which  are  tabulated  in  Table  4.1.  This  parameter  is  observed 
not  to  be  a strong  function  of  distance  along  the  body. 

Heat  transfer  and  local  Nusselt  number  results  for  Isothermal  sur- 
faces are  pi£sented  in  Fig.  4.8  for  several  values  of  wall  enthalpy 
differences  with  Mach  number  of  1.5  and  specific  heat  ratio  of  1.4. 

The  surface  is  cooled  for  values  of  wall  enthalpy  difference  less  than 
zero  and  heated  for  values  greater  than  zero.  Results  corresponding 
to  zero  heat  flux  are  for  adiabatic  surface  where  the  wall  temperature 
is  equal  to  the  free  stream  stagnation  temperature.  For  Sw  - -1.0, 
the  wall  temperature  is  at  absolute  zero  and,  thus,  there  would  be  an 
Infinite  heat  transfer  rate  to  the  surface.  Near  the  stagnation  region, 
heat  flux  is  given  by  the  function  Sq(0)  and  is  nearly  Independent  of 
distance.  A twofold  increase  of  wall  enthalpy  difference  from  1.0  to 
2.0  yields  approximately  40X  increase  in  heat  flux  for  the  stagnation 
region.  The  decrease  cf  heat  flux  with  distance  is  believed  to  be 
attributed  to  boundary  layer  separation.  Evaluation  of  Nusselt  number 
for  adiabatic  wall  condition  poses  problems  since  both  Sw  and  S'(0)  are 
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tero.  However,  results  froa  sons  preliminary  numerical  experiments 
illustrate  that  as  approaches  aero,  the  ratio  of  the  first  term  for 
S'(0)  to  Sw»  namely,  S^(0)/Sv  approaches  a value  of  about  0.76.  Thus, 
it  appears  that  the  Nucselt  number  is  defined  for  adiabatic  wall  boun- 
dary condition.  Additional  information  la  needed  to  further  define 
this  ratio. 

The  relationship  between  fluid  friction  and  heat  transfer  is 

expressed  by  the  Reynolds  analogy  parameter  which  is  shown  in  Pig.  4.9 

as  a function  of  distance  along  the  body  for  several  wall  enthalpy 

difference  values  with  Mach  number  of  1.5  and  specific  heat  ratio  of 

1.4.  Near  the  stagnation  region,  this  parameter  is  given  by  Eq.  (2.25b) 

which  includes  the  ratio  of  S /S'(0).  For  the  adiabatic  wall  condition, 

w u 

this  ratio  appears  to  acquire  a value  of  1.3.  This  then  results  in  a 
value  of  1.7  for  the  Reynolds  analogy  parameter.  This  parameter 
exhibits  values  which  are  greater  than  unity  and  which  increase  with 
wall  enthalpy  difference.  These  trends  are  similar  to  those  found  in 
other  investigations  [7,10]. 


Cf  Re*/2  Nu 
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Fig.  4.9  Reynolds  analogy  parameter  ( M<»  = 1.5,  y = 1.4) 
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5 . CONCLUSIONS 

An  analysis  has  been  developed  which  transforms  the  governing 
transport  equations  for  steady,  axisymmetric , compressible,  nonrotating 
boundary  layer  flow  about  a body  of  revolution  into  a set  of  nonlinear 
coupled  ordinary  differential  equations.  Solutions  to  the  ordinary 
differential  equations  subjected  to  specified  boundary  conditions  were 
obtained  using  a standard  numerical  integration  technique.  Results 
were  presented  for  velocity  and  temperature  profiles  within  the  boun- 
dary layer  as  well  as  skin  friction  and  heat  transfer  rates  along  the 
body. 

Several  additional  studies  are  necessary  in  order  to  completely 
establish  the  applicability  of  the  present  analysis.  First,  the 
accuracy  of  including  each  successive  term  in  the  perturbation  scheme 
needs  to  be  examined.  Neat  the  stagnation  region  of  the  body,  the 
first  term  would  be  sufficient.  However,  the  effect  of  additional  terms 
for  points  removed  from  this  area  needs  to  be  established.  Second, 
additional  results  from  this  analysis  for  other  values  of  the  parameters 
should  be  acquired  and  analyzed  for  trends.  There  is  a need  to  better 
define  the  ratio  of  S'(0)/Sw  as  Sw  approaches  zero.  Third,  comparison 
of  results  from  the  analysis  with  experimental  results  would  help  to 
establish  the  range  of  applicability  of  the  analysis.  Comparison  of 
present  results  with  numerical  solutions  of  the  governing  transport 
equations  she  ild  be  considered.  Finally,  temperatures  within  the 
boundary  layer  may  attain  sufficient  levels  where  gaseous  radiation 


can  contribute  significantly  to  surface  heat  flux.  Effects  of  radiant 
transport  within  the  boundary  layer  and  at  the  body  surface  need  to  be 
defined. 
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APPENDIX  A 


Computer  programs  employed  to  obtain  solutions  to  the  ordinary 
differential  equations  as  well  as  to  obtain  lists  and  plots  of  these 
solutions  are  furnished.  For  convenience  in  discussing  the  programs, 

Eqs.  (2.1),  (2.2),  (2.3)  and  (2.4)  are  referred  to  by  AERO,  AER1,  AER2 
and  AER3,  respectively.  The  computer  program  called  AER  used  to  solve 
for  values  of  the  unknown  derivatives  at  n ■ 0 includes  MAIN,  READ, 
RUNKUT,  FCN  and  INCON  routines.  The  purpose  of  each  routine  is  briefly 
noted  in  Table  A-l.  The  FCN  routine  is  different  for  each  AER.  Further- 
more, as  observed  by  boundary  condition  in  Eq.  (2.6),  the  routine  INCON 
is  slightly  different  for  AER3.  Initial  values  for  the  functions  as 
well  as  their  derivates  at  n ■ 0 may  be  substituted  into  the  AERL  pro- 
gram to  generate  lists  or  plots.  Routines  which  make-up  AERL  are 
briefly  described  in  Table  A-2.  Listing  of  routines  for  AERL  is  also 
supplied. 


TABLE  A-l  Computer  Program  AER 


Routine  Description 

MAIN  Controls  calling  sequence  for  other  routines. 

READ  Input  of  program  parameters  as  well  as  initial 

..lues  for  known  and  unknown  functions  at  ri  ■ 0. 

RUNKUT  Fourth-order  Runge-Kutta  integration  scheme  from 

n » 0 to  n in  steps  of  An. 
max  r 

FCN  Evaluates  functions  at  specified  value  of  n« 

Includes  perturbation  equations  for  adjusting 
guesses  [9]. 


INCON 


Adjusts  initial  values  for  unknown  derivatives 
and  checks  for  convergence. 
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MAIN 

C INTEGRATION  OF  ORDINARY  DIFFERENTIAL  EQUATIONS 
IMPLICIT  R6AL*8( A-H.U-7) 

DIMENSION  f<3,15),FO<3,15> 

C READ  INHIT  VARIABLES 

1 CALL  READ! NOD, NOR, H, ETAO, ETAM,EMAX,FO,NEQ) 
NAD J=1 

C KUNGE-XUTTA  INTEGRATION 

2 CALL  RUNKUTOOn, NOR, H, ETAO, ETAM,FO,F,NEO) 

C ADJUST  INITIAL  CONDITIONS 

CALL  I NCDN!FTAM,EMAX,FO,F,NAI)J,NEQ,.NOD) 

Wft I Tfc ( 6, 1 00 ) NAD J 
IF (NADJ.NE.OJ  GO  TO  Z 
GO  TO  l 

100  FDRMATUOX,  • INITIAL  CONDITION*  ,15  ) 

END 


READ 

SUBROUTINE  READ! NOD, NOR,H,ETAO,ETAM, fcMAX, FO, NEQ> 

C INPUT  PARAMETERS 

IMPLICIT  REAL*8( A-H.U-Z) 

DIMENSION  FO! 3,151 . T I TLS (101 

C******#*^***************)*************  *********  **********  *********** 

C NOD  NUMBER  OF  FUNCTIONS  AND  DERIVATIVES 

C NOR  NUMBER  OF  DESIRED  FUNCTIONS  AND  DERIVATIVES 

C NEO  NUMBER  OF  EQUATION  SETS 

C H INTERVAL  SIZE  FOR  ETA 

C ETAO  INITIAL  VALUE  FOR  FTA 

C FT  AM  INITIAL  VALUE  FOR  MAXIMUM  ETA 

C EM  AX  MAXIMUM  ETA 

C FO  INITIAL  BOUNDARY  CONDITIONS  AT  FTA=ETAU 

C***  ***************»**** *************** ***********  ******************* 

C READ  IN  TITLE  CARO 

REAR  ( 5, 102, E|D=99)  TITLE 
C READ  IN  ‘-'RQjR  AM  PARMETERS 

REAIM5  ,100)  NOD, NOR,  !E!J,H, ETAO, ETAM,EMAX 

c read  in  initial  values  of  f unci on s and  their  derivatives 

READ! 5,131 M (FOII , J) , J=1 ,NOD) ,1=1 ,NEO> 

C PRINT  OUT  INPUTS 

WRIT-16, 1D3)  TITLE, F0(  1,6) ,H 

WR I TE  ( 6, 1 01 ) I ( FO ( I ,J) , J=1 » NOD ) , I *1 ,NEQ) 

RETURN 

99  CAlL  EXIT 

100  FORMAT! 316,6015.7) 

101  FORMAT <3025, 16) 

102  FORMAT(IOAB) 

103  FORMAT!  Hi,l'  A8///17X,  'SW«',D15.7/ 

1 10X , * STEP  S1ZE  = '»F6.6//10X,* INITIAL  CONDITIONS*) 

END 


M O 
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RUNKUT 

SUBROUTINE  RUNKUT(NOO,NOR»H,ETAO,ETAM,FO,F,NEQ) 

C RUNGE  XUTTA  INTE3KAT I ON  SCHEME 
IMPLICIT  REAL*ft( A-H,Q-Z> 

DIMENSION  FO(3,15),F(3,15) ,FD<3, 15), FC( 3,15), AK1 (3,15), AK2I 3,15) , 
l AK 3(3,15) »AK4 (3,15) 

C ZERO  ARRAYS 

DO  7 I»1,NE0 
DO  7 J*l,NOD 
F <!,J)«0.00 
FC  ( I , J ) O»D0 
FD  1 1, J ) >0*00 
AK1( I, J)«0.D0 
AK2( t,J)»0,D0 
A<3( I, J)»0.D0 
7 A<4( I , J ) » 3*00 

C INITIAL  CONDITIONS 
16-1 

DO  1 J-l » NEO 
03  1 1*1, N30 
F( J , I ) *F3( J , I ) 

1 FC( J.I )-FO( J,I ) 

MR  I Tc ( 6, 1 01 ) 

WRITE (6, 130)  CTAO, ( F ( NEQ, 1 ) ,1-1, NOR) 

CALL  FCN(ETAO,FC,FD) 

INTE  OR AT  1 ON 

ETA-ETA3M IE-1 )*H 
00  3 J *1 « NEO 
DO  3 I = 1 » N30 
AKK  J,  I)»H*FD(  J,I) 

3 FC( J,I )-F( J, I )*0.5D0*AKl(J,I ) 

ETA*ETA*0.500*H 

CALL  FCM(€TA,FC,FO) 

DO  4 J-l , NEQ 
DO  4 1*1, NOD 
A<2(  J,  I )=HOFI)(  J,I) 

4 FC ( J.I ) = F< J,  I ) +0.5D0*AK2( J, I ) 

CALL  FCN(ETA,FC,FD) 

DO  5 J-l.NEO 
DO  5 I -1 ♦ NOD 
AK3( J, I )»H4FD( J,I ) 

5 FC ( J , I )*F(J, 1 ) ♦ AK3 ( J * I ) 

ETA-ETAO.  5t)0*H 

CALL  FCN(ETA,FC,FD) 

03  6 J * 1 , N E U 
DO  6 1-1, NOD 
A.<  4 ( J , I > *H*FD  ( J,  l ) 

6 F ( J , I ) - F ( J , I ) «-  ( AKl  ( J , I ) + 2»  DO  *(  AK2  ( J,  I )*AK3( J,I  ) )+AK4(  J,I  ) )/6.D0 
IE-IE+1 

IF(ETA.LT.ETAM)  GO  TO  2 

WR!TC( 6,100)  ETA»(F(NEQ»I)»I*1 »NOR) 

RETURN 

100  FORMAT (F7»2,5025«16) 

101  FORMAT (/3X,’ ETA' ,12X,»F» ,24X,'F*» ,23X , • F**» , 24X, »G* , 23X , » G* • ) 

END 
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PCN  FOR  AERO 


SUBROUTINE  FCN(£TA,F,FD) 

C FUNCTIONS  FOR  AERO 

IMPLICIT  REAL*0(A-H,O-Z) 
DIMENSION  F ( 3 » 15 ) , FD ( 3 *1 5 ) 
C DIFFERENTIAL  EQUATIONS 


c x 


FD  ( i 
FD(  1 
FDIl 
FO(l 
FJ(1 
PERTUR 
FDIl 
FDIl 
FDU 
l-FIl 
FDIL 
FOIi 


C V PERTURBATION 


Ff)(  1 
Fu(  * 
FOIl 
1 >-F  ( 
FDI  * 
FDIl 
RETURN 
END 


1 l«F€ It  2) 

2 )  *F ( 1 * 1) 

3) »-’.D0*F<l ,1 >*F(1,3)*F(1,2)*F(1,2)-F( 1,4)-1.D0 
4MFII.5) 

5) 3-2»D0*F(l,l)*F(l#5l 
AT10N 

6) «F< 1,7) 

7>«F{ 1,8) 

8>«-2*D0*(F(l  ♦ 6)*FU.3)*F(i»l)*F(l,8)  )*2.DO*F(l»2)*F(l,7> 
9) 

9>«F(1,10) 

13)3-2»f)0*(F(l*6)*F  J1.5>+F(1,U*F  (1,101) 


11)  *F (1,12) 

12)  »F (1,13) 

13) »-2»OO*(F(l,ll)*F(l,3)*F(l.l)*F(l,13))*2.D0*F(l,2)*F(l,12 
*14) 

14) »F(1,15) 

15) 3-2.D0*(F(l,il)*F(l,5»*F(i,il*F(i,15) ) 
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rCN  FOR  AER1 

SUBKiiUT  1 'JE  FCN(ETA,F,FO) 

IMPLICIT  REAL*fl ( A-N,  0-  Z I 
DIMENSION  F<  3,15),FO(  3,15) 

C AERO  EQUATIONS 

FDU,1)*FC1,  2) 

FOC i,?)*F ( 1, 3) 

FD(  1, 3)«-2.0:*F(l,i)*F<l  •3)+F(l,2>*F(l,2)-F(l,4)-l.D0 
FD< i,4)»F(l,5) 

F0(  i,  5 D0*FC1,1)*F  (1,5) 

C AER1  FOUAT  1 0 S 
C 01 FF  £<EN 11 4L  EQUATIONS 
FD(2,1  I a F ( 2 » 2 ) 

FD(2,2)*FI2, 3) 

FD<2,3)»-2.DD*FCl,l)*F<2»3>*4.D0*FU,2)*F<2,2l-4.D0*F(i,3)* 

1F( 2, 1 ) “2*  03*F ( 1, 1 ) *F 1 1 , 3 )-4. DO*I 1«D0+F 11, A) I -F (2,4) 
FD(2,4)*F<2,5) 

FU(2,5)»-2.00*F( 1 , 1 ) *F  (2 .5 ) ♦ 2.  DO*FU  . 2 > *F<2  »4  > -4,D0*F(  2 , 1)* 
2F(l,5)-2.0D*F(l,l»*F(l,5) 

C X-PERT UR RATION 

FD( 2, 6 ) *F ( 2*  7 ) 

FD(2,7)=FI2,B) 

F0(2,8)a-2.D:!*F(1,1)*F(2,8)+4.D0*FC1,2)*F(2,7)-4,D0*F(1,3)* 

3F(  2,6)-F( 2*9) 

FD(2,9)*F(2,10) 

FO ( 2 » 19) *-2«  00*F (1,1)*F(2,10)+2«D0*F(1*2)*F(2,9) -4. DO*F 11,5)* 
4F(  2 , 6 ) 

C Y-°ERT  ORATION 

FD( i , i 1 ) »F ( 2 , 12) 

FD(t,12)«F(2,23) 

FD(2,13)«-2,00*FI1,1 ) *F ( 2, 1 3 ) *4, DO*F 1 1 , 2 )*F( 2, 12 )~4. 00*F 1 1 ,3 ) * 
5F(  2, 11  )-Fj  2,  ’.4) 

F0( 2 , 14) *F ( 2 ,15) 

FO(2,15)--2. 00*F(l,l)*FI2,15)*2,00*FIl,2)*FI2,14)-4«00*F(l,5)* 
6F ( 2,11) 

RETURN 

END 
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FCN  FOR  AER2 

SUBROUTINE  FCNI £TA*F,FD) 
lMPLi:iT  REAL*8IA-H,0-Z) 

DIMENSION  F( 3,151 ,FD<3, 151 
C AERO  EQUATIONS 

FD(1,1)=F(1,2) 

FD(i,2)*FU,3) 

FDIi,3>=-2.D0*F(i,l>*FU#3)  + F<lt2l*FIl,2>-F(l,4l-WD0 
FD<1,4)=F<1,5) 

FD<  :,5)=-2.DC*FIltn*FU  *51 
C AER2  EOIIAT  IONS* 

C OIFFPRFNTIAL  EQUATIONS 
FOIZtl ) *F ( 2 * 2 ) 

FD(  2 t 2 ) *F  ( 2 , 3 ) 

FD(2,3>=-2.DC*Fll,l>*F(2.3>+6.D0*FIl,2)*FI2t2l-6.D0*Ftl,3)*F(2,l> 

1- 4.D0*F(1,1  >*F(l,3)-6.D0*Cl.D0+F<l,4n-FI2,4> 

F0(2,4)=F(2,5) 

F[M2»5)*-<.r)C*F<l,!>*F(2»5)M.OO*F(l,2)*F(2,4)-6.DO*F(lf  5)*F(2,1I 

2- A.D,J«F(l,:)*F(lt5) 

C X-PERTURRATIO  1 

FD(2,6)=F(2,7) 

FDI  2 ? 7 ) »F  ( 2 1 B ) 

FD(2t8)  «-2.0r*F(l,l  ) *FI2,Bl  + <>.DO*F(l,2)*F(2»7)-6.DO*F(l,3)*Ft2,6I 
7-FI2.9) 

F0( 2 » 9 ) *F ( 2 * 10) 

FDI2,IO)»-2.00*F  (1 ,1  )*F(2tlOI+A.r)0*FlI,2)*F(  2,9>-6.D0*F<  1,5)* 

3F(  2 » 6 I 

C Y-PERT UR RATION 

FDI 2 1 11 ) *F  ( 2 1 12 ) 

FD(2»12)>F(2tl3) 

FO(2,13)«-2.nO*F(l,l ) *P< 2, 13)*6,00*F < 1, 2 J*F ( 2,12 l-6«D0*F (1 » 31* 

5F( 2 » 1 1 )“F ( 2 » 1 4 ) 

FD(2,14I*F(?,15> 

Ft)  12*15)* -2  • DO*F  ( 1,  .1 ) *F(  2, 1 5 i+4.D0*F(  1 , 2 >*Ff  2,14>-6.D0*FIi*5)* 

1F( 2, 11  ) 

RETURN 

END 
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FCH  FOR  AER3 

SUBROUTINE  FCN(ETA,F,FD) 

IMPLICIT  REAL*8IA-H,0-Z> 

DIMENSION  F(  3,15)  ,FOI3,15> 

C AERO  EQUATIONS 

FDC  ,l)«Flli2) 

FDCi,2»«F(l,  3) 

F0(i,3)»-2.D0*F(1.1)*FI1 , 3 ) *F (1 , 2 ) *F ( 1 , 2 J -F 1 1 , A ) -1 .00 
F()(1,4»«F(1, S» 

F0(1 ,5)*-?,OC*FIl,l)*F(l,5> 

C AER1  EQUATIONS 

FD(L*1)*FI2»2) 

FD(?»2)*F(2,3) 

FD(2,3)»-2.3C*F( 1 , 1 ) *F < 2 , 3 »+4. DO*F t 1 , 2 > *F < 2 , 2 > -4,D0*F 1 1 , 3 ) * 
1F(2, 1)-2.D0*F(1,1I*F<1 ,3)-4.DO*ll.DO*Fll,4n-F(2,4) 
FD(2,4)*F(2»5) 

FU(2,5)*-2*DC*Hl,l)*F<2,5  5*2»D0*Fti,2  ) *F ( 2, 4) -4,D0*F ( 2 , 1 ) * 
2F(1,5I-?.DC*F(1#1)*F(1,5) 

C AER3  EQUATIONS 
C DIFFERENTIAL  EQUATIONS 
FD ( 3 » i ) 3 F ( 3 * .?  ) 

FD<3,2)«F(3,3> 

FD(3'3)--2.D:«F(l,l)*F<3t3)+6.DG*F(lt? ) *F < 3, 2 ) -6.DG*F ( l , 3 ) * 

IF  I 3,1 ) -F<  3,4  >*3.D0*F(2,2  )*F(2,2)-4.D0*F(2, i ) *F ( 2 , 3 > -2. DO*F < 1 , 3 > 
2*F(2,1  )-2.D0*F(  1,1)*F(2,3)-3.D0*(1.D0*F(1,4)  ) ♦2.00*F  (1 , 1 ) * 

3F(1 ,3)-4.D0*F<2,4) 

FDI3,A)*F(3t‘>) 

FD( 3,5 )=-2.DC*F( 1,1) *F ( 3, 5 ) «-4.D0*F ( 1 , 2 ) *F ( 3, 4 ) +2.D0*F ( 2 , 2 ) * 

1F( 2,4>-42D0*F<2tl)*F  < 2 , 5 i -2. UO*F ( 1 , 1 ) *F ( 2 , 5) -6. DO*F ( 1 , 5 ) * 

2F(3, 1 »-2.D3*F <1,5»*F<2.1 > + 2. DO*F ( l , 1 ) *F { 1 , 5 1 
C X-PERTURfATlON 

FD  l 3 , 6 I =»F  ( 3,  7 ) 

F0( 3 , 7 ) *F  ( 3,  H ) 

FD<3,B)3-2.0r*F(l,l )*F 1 3,8 >+6.D0*F 11,2) *F( 3,71-6. DO*F( 1,3)* 
1F<3,6)-F<3,9) 

FD(  3, 9 ) »F  ( 3,  ! 0) 

FD(3,10)«-2.D0*F(l, 1) *F( 3 t 10 ) ♦4.D0*F ( 1 ,2)*FI 3,9)-6.00*F( 1,51* 
1F( 3,6) 

C Y-PERTUrtHATI DM 

FD(3,11»«F(3,12) 

FD(3,12)«F(3,13) 

FD(  3,1  3U-2.DC*F  (1,  1 )*F(  3 ,1 3 ) *6. DO*F  C 1 , 2 ) *F<  3,32)-6.00*F <1 ,3)* 
1F( 3, 1 1 ) -F ( 3, ! 4) 

FD<3,14)«F<  3,’S) 

FD 1 3, 1 5) «“2#  DO*F ( 1 , 1 1 *F ( 3 , 1 5 I +4. DD*F ( 1 , 2) *F( 3, 14)-6.D0*F (1,5)* 
1F( 3, 1 1 ) 

RETURN 

END 


ijEgaasarasE 
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IHCOM  FOR  AERO,AJ5R1  AND  AER2 


SUBROUTI  NE  I NCON ( n TAM  t EMAX*FQF»FF  ,NADJ , NEQ, NOD) 

C ADJUST  INITIAL  CONDITIONS 

IMPLICIT  *EAL*8< A-H.Q-Z  > 

DIMENSION  FO<  15)  * F C 1 5 1 » FOF  ( 3 1 1 5 ) iFF  1 3, 15 ) 

03  I I-1.N3D 
FOm»FDF(Xt=Q,I » 

1 F(  I l*FF( NEQ, I ) 

All«F(7)*F(7)*F{9>*Fm*F(8)*FIimF(10)*FllO) 

A12-FI 7)»F(12)*F(9) *F( U >+F ( 8 ) *F ( 13 ) *F ( l 0 ) *F <15 ) 

A21* A1 2 

A22“F(1?)*FC12)+F(14)+F(14)*F(13)*F(I3)*F{15)*FI15) 

B1 »-  I < F I 2 > -1 . DO ) *F  ( 7 ) ♦ F U ) *F  < 9 ) *F  <3 ) *F  (8)  *F  1 5 ) *F  ( 1 0 » ) 

B2**( ( F( 2 ) -1 »D0 I ♦F ( 12)*-F|4)*F{ 14I»F(3I*F(13I*FC5)*F( 15) ) 

0EX*AI1*-A22“A12*A21 

DEX*(A22*rtl-A12*B2)/DEN 

OEYs( A11*32-A21*B1 l/UEN 

F0( 3 ) *F0 ( 3 ) ♦ DeX 

F3(5)*F3(S)  + r>fcY 

FOFINEQ, 3)«FO(3) 

F3F(*I  50,51  “F  3(5) 

c convergence  checks 

WRI TL (6, 1 30)  DEX.OEY 

IF(rjA9S(  DEX/FOI3)  ).GT.,10-12.OR.DABS(0EY/FOC5)  ).GT.. 10-12)  RETURN 
E*<F(2)-1.30)**2+F<<,)*F(A)*Fm*Fm«-F(5)*F(5) 

HR  I TE ( 6, 1 01 ) E 
ETAM-2.00*ETAM 
HR  I TE ( 6, 1021  ETAM 
IF(ETAM.Lc.EMAX)  RETURN 

etam*emax 

NADJ  »0 
' RETURN 

100  F3RMAT113X, 'CHECK  OEX  AND  OEY»,2Di5.7) 

101  FORMATIIOX, 'CHECK  E ' • 1 OX ,015, 7) 

102  FORMA) (10X, 'CHECK  ET AM • , 7X , 01 5. 7 ) 

END 


XNCON  FOR  AER3 


SUBiOUT I N£  I \'CON(ETAM,EMAX,FOF,FF,NADJ,NEQ,NOD) 

C AOJUST  INITIAL  CONDITIONS 

IMPLICIT  REAL*BIA-H,0-2> 

01  MfNS I ON  F0(15)  ,F(i;>)  ,FOF(  3, 15>,FF(  3,15) 

00  1 I «1 , NOD 
FOm-FOHNEQ.I) 

I F( I)»FFlN£Ut I I 

All. FI  7)*F<7»*FI9-*F«9**FI8)*FI8UF(10)*F(1DI 
A12-FI  7>*F(l2)4F(<n*i-(l4)4F(a)*FU3t*F(lO)*FU5> 

A21.A12 

A22»F(12)*F(12)+FI14)*FI14)4'F(l3l*FI13)*F(15)*F(i5) 

Bi.-I F ( 2 ) *F ( 7)4F( A)*FI9) ♦F(3)*FC8»»F(5)*FIIDI) 

B2»-IF(2»*Fll2)4F(A)*FllA)4F(3l*F<Ul4FI5>*F‘15n 

OEN«a11*A22-AJ2*A21 

DEX.<  A22*Bl-Ai2*B2)/t)EN 

OEY-I A11*»2-A21*B1)/DEN 

FO I 3 1 «F0( 3 ) 4DE  X 

FO 1 5 1 . F0( 5 ) ♦OE Y 

FOF ( NE  S, 3 ) «F0 ( 3 ) 

FOF(NEQ,5)«FO(5) 

C CONVERGENCE  CHECKS 

HR  I T E < 6 , 100  I OEX » OEY 

IF(0ABS(0EX/F0(3) I.GT..10-12.0R.0ABSIDEY/F0I5I UGT.. 10-121 
E«F(2)*FI2)4Fm*FU)4F!  3)  *F  I 3l*F  15 1 *F  1 5 ) 

HRI TE ( 6, 101  ) E 

E T AM*  2 • D 0 * E T AM 

HRITEI6.102)  ETAM 

IFIETAM.LE.EMAX)  RETURN 

ETAH.EMAX 

NAOJ.O 

RETURN 

100  FORMATUOX, 'CHECK  DEX  AND  DEY'f2D15,7> 

101  FORMAT  1 10<» 'CHECK  E • , 10X ,01 5. 7 ) 

102  F0RMATI10X, 'CHECK  ETAN* , 7X, 015,7) 


TABLE  A- 2 


Computer  Program  AERL 


Routine 

Description 

MAIN 

Input  of  program  parameters  and  boundary  conditions 

at  n ■ 0 for  all  functions  and  their  derivatives. 

Lists  results  and  controls  all  routines. 

RUFKUT 

Fourth-order  Runge-Kutta  integration  scheme. 

FCN 

Evaluates  functions  at  specified  n 

LIST  OF  AER 

IMPLICIT  REAL*8(A“H,0-Z) 

COMMON /DAT  A/G(5*1610) 

OINENSION  F(3,5>,FO<3,5> ,T1TLE(10) 

RE  AO  IN  INITIAL  VALUES 

RE AO! 5 * ICO* G N0«99 ) TITLE 

REA0I5* 1 D 1 ) NOR,NEQ,NPR.H,ETAO.ETAM 

READ(5,ID2)! (FOCI, J) ,J*1 *NOR) ,1*1 »NEU) 

RUNOE-KUTTA  I NTEGRATIOM 

CALL  RUNXUT!  .NOR,  N0R,H,  ET  AO,  E TAM»FOt F ,.NEQ  > 

PRINT  OUT  RESULTS 

NOATA*(ETAM-ETAO)/rt*2.001D0 
! PR*1 

WRITE (6*103)  TITLE 
002  I «1 » NDATA 
IFII.NE.IPR)  GO  TO  2 
IPR«IPR»NPR 
ETA=FTAO*H*( t-1) 

WRITE! 6,10A)  ETA, (S( J, I ) ,J»1,N0R) 

CONTINUE 
GO  TO  \ 

9 CALL  EXIT 

00  FORMAT ( 1DA8) 

01  F0RMAT(3IS,3D15.7) 

02  FORMAT!  3025.16)  ....... 

03  FORMAT!  HI,  10A8/3X,  ‘ETA*  ,8X,  »F'  ,1*X,  'F*»  ,11X, ' F**'  ,13X,  *G*  tlAX, 

1 *G*  • ) 

0*  F0RMATIF7. 2,5015. 7) 

END 


M ft  0»  *-»  O 
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RUNKUT 

SUBROUTINE  RUNKUT < NOD, NOR  ,H , ETAO, ET AH, FO,F,MEQ) 

C RUNCE  KUTTA  INTEGRATION  SCHEME 
IMPLICIT  REAL*8(A-H,0-Z) 

COMMOM/DATA/GI 5,1610) 

DIMENSION  POO,  15)  ,F ( 3,1 5) ,F0<3, 151 ,FC( 3,15) , AK1 (3,15 1, AK21 3, 15) , 
1 AK3( 3,15) ,AK4(3,15) 

C ZERO  ARRAVS 

DO  7 1*1, NEQ 
DO  7 J*l, NOD 
F ( I , J ) *0»D0 
FC  (l,J)*0.00 
FD  <I,J>*0.00 
AKKI,  J)*O.DO 
AK2(  I, J)*0.00 
AK3< I » J)*0»00 
AK4( I, J)*0.00 
INITIAL  CONDITIONS 
IE*1 

DO  i J *1 » N£Q 
DO  1 1*1, NOD 
F(  J , I ) *F0( J,  I ) 

FC ( J, I )«F0( J,I  ) 

DO  8 1*1, NOR 
G( 1 , IE ) *F ( NEQ, I t 
CALL  FCN(ETAO,FC,FD) 

INTEGRATION 

ETA*STA3M IE-1 )*H 
03  3 J *1 , NEQ 
DO  3 1*1, NOD 
AK1I J,I»*M*FD( J,I) 

3 FC( J,I  >*F< J, I H-0.5D0*AKl(J,I ) 

ETA=ETA*0.50D*H 

CALL  FCN(ETA,FC,FD) 

DO  A J * 1 , N c Q 
DO  4 1*1, NOD 
AKZI J, I J=H*FD( J,I ) 

4 FC  I J,  I )*F(J, I) ♦0*5D0*AK2 ( J, I ) 

CALL  FCNI ETA , FC, FO ) 

DO  5 J*1,.N?Q 
DO  5 1*1, NOD 
AK3 ( J, I )»H*FO( J,I  ) 

5 FC ( J , I ) *F ( J * I)+AK3(J,l) 

ETA*ETA  + 0.5D'MH 

CALL  FCN(ETA,FC,FD) 

DO  6 J*l" , NEO 
DO  6 I * 1 , N )0 
AK4U,  I )«H*Fn(  J,I) 

6 F(  J,I  )*F(  J,I » ♦ | AK1  ( J,  I ) + 2.D0*(  AK2(J,1  >«-AK3(  J,I>)+AK4f  J,  I ) )/6*D0 
IE*  IE* 1 

DO  9 1*1, NOR 
9 G( 1 , IE ) “F ( NEQ, I ) 

IF(ETA,LT.ETAM)  GO  TO  2 

RETURN 

END 
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FCN  FOR  AERO 


SUBROUTINE  FCNI £ TA»  F ,FD) 

C FUNCTIONS  FOR  AERO 

I I C I T REAL*«<  A-H,U-Z  ) 

OIH£5.'4SIOM  F ( 3, 15), FD (3,15) 

FIX*,  U»FC1,  21 
FDU,2 1.FIl,  3i 

FO(i,3i«-i.OOtF(  1,1>*F(1 ,3H-F(lt2)*Fa«2)-F(l>4l-UO0 
FOCi,A)«F(l,S) 

FOIi , 5)»-2.DO*F< 1,1 )*F  U ,5) 

RETURN 

END 


FCN  FOR  AER1 

SUBROUTINE  FCN(£T4,F,FD) 
implicit  REAL*8(A-H,0-Z> 

DIMENSION  F( 3,15) ,FD(3,15) 

C AERO  EQUATIONS 

FDI i » 1 ) *F I 1 , ? ) 

FO( 1 f 2 ) -F ( 1 , 3 ) 

Fo!j!J!rF^i??r<1,u*F<l,3,+Hi,2,*,i<i,2,"F,i’4,"i/*Do 

FOIL, 5)*-2.DC*F<  1,1  I*FU,S) 

C AER1  EQUATIONS 

FDI2 , 1 )*F  < 2 1 2 ) 

FO(2v2)*F(2,3) 

* 1,1  ,*F(l*3,“4*U0*<i*°0+F(l»An-FI2,A) 

RETURN 

END 
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FCN  FOR  AER2 


SUBROUTINE  FCN  ( ETA ,F  » FD) 

IMPLICIT  REAL*8(A-H,0-Z> 

DIMENSION  F(3,15),FD(3,15) 

C AERO  EQUATIONS 

FD<1,1 »«Ftl,Z) 

FDtl.2)»Fti, 3) 

FD(1,3)--2.D0*F(1,1)*F<1 ,3 1 + F 1 1 , 2 > *F ( 1 , 2 I -FU , 41-1. DO 
FD  ( i * 4 1 «F  t 1,5) 

FDt 1 , 5 ) *-2. D0*F  ( 1 , 1 I*F  (1  ,5 1 
C AER2  EQUATIONS. 

FO<2,l)*F(2,?> 

FD<2,2>*F(2,3> 

FDt  2, 3 I *-2»DO*F  tl , 1 )*F (2 .3 l+6.D0*Ftl,2 1 *F 12*2 l-6*00*FU t 3)*F! 2, 1 1 

1- 4«OP*Ftl,l)*F<l,3>-6«DO*tl«DO*F<lf4)  I -Ft  2*4) 

FD<2,4)*F(2,5) 

FDt2,5)»-2»D0*Ftl»l)*F(2,5)+4*D0<-F(l,2>*Ft2f4)-6«D0*Fti,  5>*Ft2,l) 

2- 4. DC* F tltl)*F(l»5l 
RETURN 

END 


FCN  FOR  AER3 

SUBROUTINE  FCN t ETA ,F , FO) 

IMPLICIT  R£AL*8(A-H,0-Z) 

DIMENSION  F(3fl5)»F0(3tlSI 
C AERO  EQUATIONS 

FDUtl  )*FU.i  ) 

FD ( 1,21 *F( 1,3) 

FO(l,3)«-2.DO*F(l,l)*Faf3I  + F<l,2)*F<l,2l-Ftl,4)-l.DO 
FO ( 1 , 4 ) *F  ( 1,5) 

FD(i,5i*-2.D0*F(l,n*F(l  ,5) 

C AER1  EQUATIONS 

FD  ( 2, 1 ) =F  < 2,  2 ) 

FD  ( 2 , 2 ) =F  ( 2 , 3 ) 

FDt2,3)«-2.0n<  F(l,l ) *F ( 2 , 3 I +4.D0*F 1 1 , 2 ) *F ( 2 , 2 I -4. D0*F 1 1 , 3) * 

1F(  2,1 1-2.03*1  ( l,l)*F(lt3)-4.DO*(UDO+F(l,4))-F(2,4) 

FD(2,4)*F( 2,5) 

FO<2,5  ) (1,1)  * F (2  , 5 ) ♦ 2# DO*F  1 1 » 2 )*F(2,4)  —4  »D0*F(2,  1)* 

2F< 1,5)-2.D3*F< i * 1 ) *F ( l ,5 ) 

C AER3  EQUATIONS 

FD<3,1 )=F(3,2) 

Fl)(3,2)*F(  3,3) 

FDt  3,3)«-?.DG*t :(1,  1 ) *F  ( 3,3  >♦  6.00*  F 1 1 , l*F(  3,  2 ) -6. I>0*F  ( 1 , 3 ) * 

1F<  3,1)-F(3,4)«  3.D0-  f (2,2)*F(2  , 2 I - 4.  OC<'F  ( 2 , 1 1 *F  I 2 , 3 >-2.  DO  *F(  1 , 3 ) 
2*M2,1  )-2.D0*F  (1,1  )*F(2,3I-3.D0*(  l.DO  + FU,  4)  )+2.D0*F(  1,1  )* 

3F( 1,31-4. DC*F(2, 4) 

FD(  3,4 )«F( 3, 5) 

FI>t  3,5)«-2.Df -N  (1,  J)*FI3,5  )*4.D0<  Ft  1,2  ) *F  ( 3,  4 ) *2.  DO*F  ( 2.  2 ) * 
iF(2,4)-4.I»0*F(/,l>  *F  (2,5  )-2,00*F(  1,1  KF|  2 , 5 1 -6.  DO*F  ( 1 , 5)  * 
2F(3,l)-2.r)0*F(l,5)'-F(2,l)*2.m*F(l,ll<F(l,S) 

RETURN 

END 


